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^\ , Local density of states (LDOS) in the triangular vortex lattice is investigated based on the quasi- 

^\ ' classical Eilenberger theory. We consider the case of an isotropic s-wave superconductor with the 

t-H | material parameter appropriate to NbSe2. At a weak magnetic field, the spatial variation of the 

LDOS shows cylindrical structure around a vortex core. On the other hand, at a high field where 
the core regions substantially overlap each other, the LDOS is sixfold star-shaped structure due to 
the vortex lattice effect. The orientation of the star coincides with the experimental data of the 
scanning tunneling microscopy. That is, the ray of the star extends toward the nearest-neighbor 
t-H ' (next nearest-neighbor) vortex direction at higher (lower) energy. 
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PACS numbers: 74.60.Ec, 74.25.Jb, 61.16.Ch 
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1 I. INTRODUCTION 

5h ■ 

Since the success of observing the vortex core image by Hess et al. in the scanning tunneling microscopy (STM) 
experiments, it has been possible to investigate experimentally the local density of states (LDOS) around a vortex core 
in type II superconductors UU This investigation has revealed a rich internal electronic structure associated with a 
vortex core, and should be important for understanding the conventional and unconventional superconductors. Then, 
to extract further information from the experimental data, the detailed calculation of the LDOS is expected from the 
13 . theoretical side. 

Hess et al. have done a series of beautiful STM experiments on a layered hexagonajLcompound 2H-NbSe2 (T C =7.3K) 
to reveal the detailed spatially resolved electronic structure around a vortex core.EJia They provided direct images of 

i^i - individual vortices and the flux line lattice. Their experiment is the first one detecting the.-c|uasiparticle state bound 
to the vortex core, which was theoretically predicted by n Caroli, de Gennes and Matricon.Q Further, the remarkable 
. result is the sixfold star-shaped LDOS around a vortexoO Their results are summarized as follows when the magnetic 
field is applied perpendicular to the hexagonal plane: (1) The LDOS for quasiparticle excitations has a sixfold star- 
shape centered at a core. (2) The orientation of this star depends on the bias energy. At zero bias, the "ray" of the 
star extends away from the nearest-neighbor direction where the conventional 60° Abrikosov vortex lattice is formed. 
Upon increasing the bias voltage, the star rotates by 30° and the ray extends to the nearest-neighbor direction. (See 

_ ' Fig. 4 in Ref. || or Fig. 1 in Ref. [|.) (3) In the intermediate bias voltage, a ray splits into a pair of nearly parallel 

\& [ rays, keeping its direction fixed (see the STM image for 0.48 mV in Fig. 1 of Ref. |j). 

ON . Recently, on the other hand, Maggio- Apple et al. have succeeded in observing the STM image of the vortex cores 
on a high-T c superconductor YBa2Cu307.@ One of the point in analyzing their STM image of vortex is how the gap 
anisotropy (d x 2_ y 2-w&Ye in the high T c superconductor) affects the LDOS around a vortex. 

To understand these experimental results, the concrete form of the LDOS structure is expected to.be calculated 
1 , [ fronipthe theoretical side. So far, the LDOS around a vortex was calculated by Gygi and SchluterJj'Ej and Shore 
et alS^ from the Bogohuhoy-de Gennes (BdG) equation. The calculations of the LDOS based on the quasi-classical 
Eilenberger (QCE) theoryEI were done by Kleinjlj and Ullah et al&3 for an isotropic s-wave superconductors. While 
these calculations showed that the vortex image observed by STM is due to the quasiparticle state bound to the vortex 
core, they investigated only the case of an isolated single vortex, and based on the assumption of the cylindrically 
• i-h , symmetric vortex core structure, rifl-4 d x 2_ y 2-wetve superconductor, the LDOS around an isolated single vortex was 
calculated by Schopohl and Maki,Ej'E£l and Ichioka et alxA However, the LDOS in the vortex lattice case has not been 
calculated so far. 

On the other hand, the STM experiments associated with the vortex core are usually performed at high magnetic 
fields, where the distance between vortices is short and the overlap of the vortex core with that of the nearest- 
neighbor vortices can not be neglected. In this situation, the LDOS around a vortex core is expected to break 
cylindrical symmetry and show sixfold symmetric structure when the vortex lattice forms an triangular lattice. The 
purpose of this paper is to calculate the LDOS in the triangular vortex lattice by using the QCE theory, and clarify its 
sixfold symmetric structure. It is expected to be observed at a high magnetic field. For the case cylindrical symmetry 
is broken, the QCE approach is more suitable than the BdG approach. 
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As for the sixfold symmetric structure of the LDOS observed by Hess et aZ.,i~i Gygi and SchluterSH discussed 
it on the basis of the BdG theory by introducing a sixfold symmetric perturbation term and using their numerical 
solutions of the cylindrical symmetric case. While they explained the above-mentioned experimental features (1) and 
(2), the origin of their perturbation term is not clear, and it is uncertain whether the sixfold symmetric term can 
be treated by a perturbation theory. We note that their theory cannot determine the absolute orientation of the 
star relative to the vortex lattice configuration, since in their theory the orientation is determined by the sign of the 
perturbation term, which is given as an assumption. 

As for the origin of the sixfold symmetric vortex structure, the following three possibilities are enumerated within 
the weak coupling BCS theory; the effect of the vortex lattice (that is, the effect of nearest- neighbor vortices), the 
effect of a sixfold symmetrically anisotropic s-wave pairing, and the effect of the sixfold symmetrically anisotropic 
density of states at the Fermi surface. In this paper, we focus on the effect of the vortex lattice by calculating the 
LDOS naively in the vortex lattice case. The contribution of this effect should be clarified before considering other 
anisotropic effects. The contribution of the other anisotropic effects is investigated elsewhere.Ej 



The quasi-classical calculation in t 
II superconductors such as Ta or Nb. 



vortex lattice case was so far performed by Klein in the case of 1ow-k type 
He calculated the spatial variation of the pair potential and magnetic field, 
where the Eilenberger equation is solved self-consistently in the Matsubara frequency and he succeeded in solving 
it by eithep-the so-called symmetry method or the so-called explosion method.ElO He also tried the calculation of 
the LD0S,E3I23 where the Eilenberger equation is solved in real energy instead of the Matsubara frequency by using 
the self-consistently obtained pair potential and vector potential. In this real energy case, he could calculate only 
the momentum-resolved LDOS for specific k (the relative momentum of the Cooper pair) direction because "A the 
symmetry method was used. In this paper, we calculate the LDOS following the method suggested by Klein.Ejllj'Ej 
However, as we succeed in solving the Eilenberger equation in the real energy case for arbitrary k directions by using 
the explosion method, we can calculate the momentum-resolved LDOS for arbitrary k directions. Then we can obtain 
the LDOS integrated over all k directions. 

The LDOS in the vortex lattice is an important physical quantity, since it can be observed directly by the STM 
experiments. And further, it can be a clue of estimating the transfer of the quasi-particle bound— states between 
vortex cores. This transfer leads to the band structure of the bound states as suggested by Canel.Bil If the enough 
transfer exists, it makes the de Haas- van Aiahen (dHvA) oscillation possible even in the superconducting state, which 
is observed in the materials such as NbSe2c3 and YBa2Cu307.c3 

In our|-|Calculations, we consider the same situation as that Hess et al. performed their STM experiments on 
NbSe2 B~U We use the material parameter appropriate to NbSe2 , and assume that the Fermi surface is two dimensional, 
which is appropriate to NbSe2- The magnetic field is applied perpendicular to the hexagonal plane, that is, along 
the z axis. The case of the triangular vortex lattice in an isotropic s-wave superconductor is considered in the 
clean limit. The Fermi surface and the energy gap of the superconductivity are assumed to be isotropic in order to 
exclude other origins of anisotropy and to clarify the vortex lattice effect on the LDOS around a vortex. Throughout 
the paper, energies and lengths are measured in units of the uniform gap A at T = and the coherence length 
Co = v-p / Ag = 7t£bcs (i>f is the Fermi velocity and £bcs BCS coherence length), respectively. The magnetic field and 
the vector potential are, respectively, scaled by 0o/£ 2 and <fia/£, where <j>o is the flux quantum. 

The rest of this paper is organized as follows. In Sec. [n[ we describe the method of calculation based on the QCE 
theory. Section III contains num erically obtained results about the LDOS of the vortex lattice. The summary and 



discussions are given in Sec. IV 



II. QUASI-CLASSICAL EILENBERGER THEORY 

The calculation based on the Eilenberger theory is performed by following the method suggested by Klein. 
First, we obtain the pair potential and vector potential self-consistently by solving the Eilenberger equation in the 
Matsubara frequency. Next, using them, we calculate the LDOS by solving the Eilenberger equation in the real energy 
instead of the Matsubara frequency. 

In our calculation, the unit vectors of the vortex lattice are given by ri = (a x , 0), Y2 = {C,a Xl a y ). As we consider a 
triangular lattice, we set a y /a x — v / 3/2 and £ = 1/2. The microscopic magnetic field H(r) = (0, 0,H(r)) is divided 
into an external field H = (0, 0, H) and an internal field h(r) = (0, 0, h(r)), 

H(r) = V x A(r) = H + h(r), (2.1) 

where the spatial average of h(r) vanishes. Therefore, the vector potential A(r) is also divided into two parts, 

A(r) = ifl x r + a(r) (2.2) 
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in the symmetric gauge. From Eqs. (^J) and ( |2.2| ), we obtain h(r) = V x a(r). 

We consider the quasi-classical Green functions g(u) n , 9, r), f(ui n , 9, r) and /'(w n , 9, r) with the Matsubara frequency 
u> n = (2n + 1)tcT, where r is the center of mass coordinate of a Cooper pair. The direction of the relative momentum 
of the Cooper pair, k = k/|k|, is denoted by an angle 9 measured from the x axis in the hexagonal plane. For the 
quasi-classical Green functions, the Eilenberger equation is given as 



{w„ + i (d\\ + ^ A \\) }/K- 6, r) = A(r)gK, 9, r), 



(2.3) 



{»n-l(t 
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}fi(u n ,6,r) = A*(r) S K,0,r), 



(2.4) 



a {l g(u n , 9, r) = A*(r)/K, 9 7 r) - A(r)/t( u , 0, r), 



(2.5) 



;K,0,r) = (l-/(wn,fl, r)/t(w n ,fl,r)) * Re 5 (u; n , 0, r) > 



(2.6) 



where 9| | =d/dru and Am = k • A = — \Hr±_ + k • a. Here, we have taken the coordinate system: u = cos0x + sin6*y, 
v = — sin#x + cos#y, thus a point r = + yy is denoted as r = rnu + r±v. The first-order differential equations 



(2.3)-(OI) are solved along the trajectory where r± is held constant. 



The self-consistent conditions for the pair potential A(r) and the vector potential are, respectively, given as 



A(r) = VN 2ttT £ / * ^-f(ui n , 9, r), 



(2.7) 



Vx Vx A(r)-Vx Vxa(r) = 42tf V / — ~<?K,0,r), 



(2.8) 



where Nq is the density of states at the Fermi surface, V the pairing interaction, and k = (7C(3)/72) 1 / 2 (Ao/fc_BT c )fto 
0.603^0 with Ricman's zeta function C(3). The GL parameter kq is given by 



7ttC(3) 



(hv F ) 4 
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In our calculation, we use the relation 



VN 



ln^ + 2vrT ]T 



(2.9) 



(2.10) 



0<U1„<U1 C 



and set the energy cutoff lu c = 20T C . When A(r) is calculated from the differential equation (2.8), the Fourier 
transformation is employed by following the method of Klein.Ej 

We calculate the r.h.s. of Eqs. ( |2 . 7| ) and (2.8) using the quasi-classical Green functions obtained by Eqs. ( p.3[ )-( ^6[ ), 
and obtain the new v alue f or A (r) and A(r). Using the renewed pair potential and vector potential, we solve the 
Eilenberger equation (2.3)- (|2.6| ) again. Starting from the initial form, 



A(r) 



1/4 



I a y I y + y 
exp < — 7r— h p 

a x \ a y 
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xy 



(2.11) 



and a(r) = 0, we repeat this simple itera tion procedure more than 20 times, and obtain a sufficiently self-consistent 
solution for A(r) and A(r). In Eq. ( [2.11 ), the r.h.s. is the Abrikosov solution of the vortex lattice, where we use the 
relation Ha x a y /4>o = 1- The factor exp(— iirxy / a x a y ) is due to the gauge transformation from the Landau gauge to 
the symmetric gauge. We set r = (xoiVo) = ~\i v i + r 2) s ° that one of the vortex centers locates at the origin of 
the coordinate. 

On the other hand, the LDOS for energy E is given by 
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f'^ 17 dO f^ 71 d6 

N(E,r) = —N{E,v,9) = —Reg(iu n ^E + ir),6,r), 
Jo 2lT Jo 27r 



(2.12) 



where 77 is a positive infinitesimal constant , N ( E,r, 9) the angular-resolved LDOS. Typically, we choose ij = 0.01. To 
obtain g{iw n — > E + irj, 9, r), we solve Eqs.(2.3)-( |2.6| ) for r\ — iE instead of uj n using the self-consistently obtained pair 
potential and vector potential. 

In our calculations, we discretize r = ur\ + W2 in a unit cell (—0.5 < u < 0.5, —0.5 < v < 0.5) for A(r) and a(r). 
There, u and v are discretized into the 41 mesh points, respectively. When we solve Eqs.( |2.;j| )-(2.G), we need to know 
A(r) and a(r) for arbitrary r. It is given by interpolation of the value on the mesh point. It is enough that A(r) and 
a(r) are calculated within a unit cell. Once we obtain A(r) and a(r) in a unit cell, we can know A(r) and a(r) in 
other unit cells by the following relation of the lattice translation R = mr\ + nr 2 (to, n: integer), respectively, 



A(r + R) = A(r)e' ix(r ^ R) , a(r + R) = a(r), 



where 



x(r,R) 



2tt 



A(R) 



2tt 

nmn + — (H x r ) • R 

0o 



(2.13) 



(2.14) 



in the symmetric gauge. 

Since we need a lot of CPU time for the integrating process of Eqs. (2.3)-(2.5) in the calculation of the quasi- 
classical Green functions, we want tq-ahorten this process as much as possible. It becomes possible by the symmetry 
consideration, as suggested by Klein.E] Without the integrating process of all 9 and r cases, we can obtain the Green 
functions for all 9 and r. 

When one of the vortex center locates at the origin of the coordinate (ro = — 5 (1*1 + r 2)), the pair potential has the 
relation 



A(r) = -A(-r). 



(2.15) 



Considering the transformation r 



-r, we obtain the following relation from Eqs. ( 2^3 )-( 2^6 ) and ( 2.15 ), 

/K,0,r) = -/t>;,0,-r), 
/tK,0,r) = -/*K,fl,-r), 

sK,0,r) = <fK,0,-r). (2.16) 

In the calculation of the Matsubara frequ e ncy uj n or the case E = 0, once we calculate the Green functions for r in 
a half area of a unit cell from Eqs. ( |2 . 3[ ) -(2.6 ) , we can know the Green functions for all r in the unit cell from the 
relation ( 2.16 ). 

When we consider the reflection at the x axis St = (x, — y), our definition of A(r) gives the relation 



A(r) = -A* {St). 

Then we obtain the relation 

f(ui n ,e,T) = -f*(w n ,-e,ST), 

/tK,^r)=-/t> n ,-0,Sr), 
g(u n ,6,r) = g*(uj n ,-9,Sr) 

from Eqs. ([2.3|)-( |2~6| ). On the other hand, there is a relation 

A(r) = A{R n r)e mn/3 , 



(2.17) 



(2.18) 



(2.19) 



when we consider the rotation nw/3 (n: integer) around the origin of the coordinate; R n r = (xcos(nir/3) 
y sin(n7r/3), x sin(n7r/3) + ycos(n7r/3)). Then we obtain the relation 



/K, 9, r) = f{uj n ,9 + mr/3, R n v)e m ^ 3 , 
f\uj n ,6, r) = /t( Wn , 9 + mr/3, R n v)e~™/\ 
g{w n , 9, r) = g(oj n , 6 + mr/3, R n r) 



(2.20) 



from Eqs. ( J2 . 3| ) - ( J2 . 6| ) . Once we calculate the Gree n fun ctions for < 9 < ir/6 from Eqs. (2.3 )-( 2.6 ), we can know the 
Green functions for all 9 from the relations ( |2.18 ) and (|2.20| ). 

There are two methods to solve the Eilenberger equation ( 2.3 )- ( 2.6 ) ; the explosion method and the symmetry 
method. We describe them in the following. 
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A. explosion method 



We use the so-called explosion mcthodEJcIl in this paper to obtain the quasi-classical Green function 



9 if 
-9 



(2.21) 



from the Eilenberger equation ( p.3[ )-( 2.6 ) . In addition to a physical solution <7 p h, Eqs. ([2.3|)-(2.6) have two unphysical 
solutions g+ and <?_. The solutions g± explode (increase exponentially) in the directions ±k and decrease in the 
opposite directions. Even when we use the physical solution as an initial val ue, t h e un physical solutions always mix 
and become dominant during the process of the numerical integration of Eqs. (2.3)-(2.5) along a long path, .-We obtain 
g± by integrating from rii =F r A to ru, where ta(> 0) is large so that explosion takes place. It is knownEJ'Eil that the 
physical solution is obtained from the commutator of the two unphysical solutions, 



9 P h = c[g + ,gJ\, 



(2.22) 



where c is a constant determined from Eq. (2.6). In this method, we can obtain the solution for arbitrary k directions. 



B. symmetry method 



In the case of solving the Eilenberger equation in the real energy E, Klein could not succeed in the .explosion 
method. Then, he calculated the angular-resolved LDOS by means of the so-called symmetry method£3€3 The 
lattice translation R = mvx + nv^ (m, n: integer) of the quasi-classical Green function is given as 



/K,0,r + R) = /K,0,r)e^ R >, 
/t( Wn ,e,r + R)=/t(o;„,e,r)e-^ R ), 
g(ui n , 9, r + R) = g(ui n , 6, r), 



(2.23) 



where x(r,R) is defined in Eq. ( 2.14 ). When we solve Eqs. ( |2.3| )-(2.6) f° r the k direction parallel to mri + nr2, 
t he so luti on have to satisfy Eq. ( 2.22 ) as the boundary condition. Thus, we integrate the first differential equ ation s 
( |2.3D -(2.5) along the trajectory from r to r + R, and search the solution which satisfies the boundary condition ( 2. 23] ) 
by the so-called shooting method or other method. 

Compared with the explosion method, the integral path can be short in the symmetry method, especially for the 
calculation of the real energy E. Then the CPU time of the numerical calculation may be short. However, in the 
symmetry method, we can obtain the solution only for the specific k di rection parallel to mx\ + nvi- Since we must 
know the solution for all k directions to calculate the LDOS in Eq. (2.12), we cannot obtain the LDOS in this method. 



III. PAIR POTENTIAL, MAGNETIC FIELD AND LDOS 



The pair potential, the vector potential and the LDOS are calculated for the material parameter appropriate to 
NbSe2, i.e., the BCS coherence length 77 A and the BCS penetration depth 690 A.0 As an example of a low magnetic 
field, we consider the case H= 0.1 Tesla (a x = 6.4£). And as an example of a high magnetic field, we consider the 
case H= 1 Tesla (a x = 2.0£). As for the example of further high field, the cases H= 2 Tesla (a x — 1.4£) and 4 
Tesla (a x — 1.0£) are also considered. Figure [l] presents the configuration of the vortex lattice presented in figures 
in this paper. In the figures, there are 7 vortex centers, and one of the vortex centers locates at the center of each 
figure. When we consider the profile of the spatial variation, we present it along the trajectories of line OA (0° radial 
direction), line OB (15° radial direction), line OC (30° radial direction) and line AC (boundary). 



A. Pair potential and magnetic field 



The amplitude of the self-consistently obtained pair potential |A(r)| is shown in the contour plot of Fig. Q. For 
#=0.1 Tesla in Fig. | (a), the core region localizes in the small region around the vortex center in the unit cell. For 
H=l Tesla in Fig. | (b), the core region occupies the large part of the unit cell. Then the core regions substantially 
overlap each other. As seen from the contour lines, while the inner region of the vortex core (the region where 
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|A(r)| < 0.8) has cylindrically symmetric structure, the outer region (|A(r)| > 0.8) shows the sixfold symmetric 
structure. There, the amplitude along the nearest-neighbor vortex direction (0° direction) is a little suppressed 
compared with that along the next nearest-neighbor vortex direction (30° direction). In this high field case, the 
maximum of the amplitude is suppressed down to 0.95Ao. The sixfold symmetric structure of |A(r)| appears more 
clearly in the case H=2 Tesla in Fig. || (c), where the maximum of the amplitude is suppressed down to 0.83Ao- In the 
case H=A Tesla, the pair potential has similar sixfold symmetric spatial distribution, but its amplitude is suppressed 
down to 0.38Ao- 

The spatial variation of the magnetic field is obtained from the self-consistently obtained vector potential. The 
microscopic magnetic field H(r) is shown in Fig. For H =0.1 Tesla in Fig. || (a), H(r) has a cylindrically symmetric 
sharp peak at each vortex center. The variation of H(r) is in the range 0.75 < H(r)/H < 2.95. For H=l Tesla in 
Fig. U (b), this peak becomes lower and broader (0.981 < H(r)/H < 1.096). It has almost cylindrically symmetric 
but slightly sixfold symmetric structure. There H(r) extends slightly to the 30° direction in the core region. The 
sixfold symmetric distribution of H(r) is seen clearly in the case H=2 Tesla in Fig. || (c), where the peak of H (r) 
becomes further low and broad structure (0.994 < H{r)/H < 1.025). 



B. Sixfold symmetric LDOS 



The LDOS in Eq. ( |2.12| ) is calculated by using the self-consistently obtained pair potential and vector potential 
presented in Figs. and I 



As for the angular-resolved LDOS N(E,r,9) in Eq. ( 2.12 ), the results of Klein in the symmetry methodE3 are 



qualitatively reproduced by our calculation of the explosion method. In the single vortex case, N(E, r, 9) consists of 
a straight sharp ridge line parallel to the k-direction, where its distance from the vortex center increases on raising 
E. At a high magnetic field of the vortex lattice case, the ridge of N(E, r, 9) extends to the neighboring unit cells 
(see, for example, Fig. 5 in Ref. ^2|). Especially _Lt is noted that, when k points toward the nearest-neighbor vortex 
direction (9 = 0° or kio in the notation of KleinE3) , the sharp ridge in the single vortex case splits into two parallel 
ridges. There, N(E, r, 9) distributes broadly in the region between them (see Fig. 10 in Ref. |22] ). Also when k points 
toward the next nearest-neighbor vortex direction (9 = 30° or kn), N(E,r,9) slightly splits into two ridges. The 
above-mentioned splitting of the ridge occurs more clearly for higher magnetic field or larger E. At a weak magnetic 
field, the ridge of N(E, r, 9) localizes within a unit cell, and the splitting of the ridge does not occur even for 9 = 
(see Fig. 1 in Ref. E|. 



Integrating N(E, r, 9) over all 9, we obtain the LDOS N(E, r) in Eq. ( 2.12 ). The spatial variation of the numerically 
obtained LDOS is shown in Figs. Hand || for the weak magnetic field case H = 0.1 Tesla. Figure ^ shows the contour 
plot of N(E,r) for each E. Figure^ shows the profiles of N(E,r) along the lines OA, OB, OC and AC of Fig. [jj As 
shown in Fig. [l] (a), the LDOS shows cylindrical structure in the low energy case (E < 0.8), where the ridge of the 
LDOS forms a ring around each vortex center. It is the similar structure as that obtained in the single vortex case 
for s-wave pairing (see, for example, Fig. 1 in Ref. ^|). However, in the higher energy case shown in Figs. ^ (b) and 
(c), the radius of the ring approaches the boundary of a unit cell in the vortex lattice, and the LDOS becomes sixfold 
symmetric structure. 

The LDOS for the high field case H — 1 Tesla is shown in Figs. || and |?]. It presents the sixfold symmetric 
star-shaped structure even in the low energy case. At E = 0, the ray of the star extends to the next nearest-neighbor 
directions (the 30° n directions) as presented in Fig. [6] (a). The width of the ray is rather broad compared with the 
experimental dataBij At higher energy E ~ 0.6, on the other hand, the ray extends to the nearest-neighbor directions 
(the 0° directions) as shown in Fig r M.(c). The orientation of the star relative to the vortex lattice is consistent with 
the experimental data by Hess et alnu Therefore, the experimental features (1) and (2) mentioned in Sec. |, including 
the orientation of the star, are qualitatively reproduced without further assumptions. As presented in Fig. ^ (b), at 
the intermediate energy of the 30° rotation of the star, the amplitude of the ray in the 30° directions decreases and 
the new ray starts to extend in the 0° directions with increasing E. This behavior is in discord with the experimental 
data. In the experiment on NbSe2, the LDOS shows the split parallel ray structure at the intermediate energy (the 
feature (3) mentioned in Sec. Q). 

The experimental data around E ~ 1 can be qualitatively reproduced by our calculation. For E ~ 1, there are little 
LDOS distribution around the vortex core region, and the LDOS distributes around the boundary of the Wigner-Seitz 
cell of the vortex lattice. Therefore, the vortex core is detected as a dark object in the STM image. Even in the case 
the vortex core localizes in the narrow region around the vortex center, the dark STM image of the vortex core for 
E ~ 1 has the size comparable to the unit cell, as shown in Figs. ^ (b) and (c). In the LDOS distribution around 
E ~ 1, we pay our attention to the value at the boundary of the Wigner-Seitz cell, that is, on the line AC in Fig. 
[I]. For E < 1, N(E : r) is large at the point A compared with the point C, as seen from Fig. |^ (d). On the other 
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hand, the peak of N(E,r) shifts to the point C on raising E as seen from Fig. || (e), which is consistent with the 
experimental STM image for 1.2 mV in Fig. 1 of Ref.H 

The spatial variation of the LDOS structure in Fig. gcan be explained as follows. Without the vortex lattice effect, 
the LDOS forms the cylindrical structure which is shown schematically as white rings in Fig. |[ It has the ridge on 
a ring around each vortex center, and the LDOS distributes at the outer side of the rings. To this distribution, the 
vortex lattice effect suppresses the LDOS along the common tangent lines on these rings (lines in Fig. ||). These lines 
runs toward the 0° and its equivalent directions (the nearest-neighbor vortex directions). This suppression is due 
to the splitting of the ridges in the angle-resolved LDOS N(E, r, 9) for 9 of the nearest-neighbor directions. In the 
contribution to the LDOS, the ridge of N(E, r, 9) along the lines of Fig. || becomes broad and low-height distribution 
due to the splitting (see Fig. 10 in Ref. g2|). Therefore, it is seen that the LDOS in Fig. || has small value at the 
points along the lines. The small suppression is also shown along the tangent lines running toward the 30° and its 
equivalent directions (the next nearest-neighbor vortex directions). Considering the fact that the radius of the ring 
in Fig. ^ increases on elevating E, we can understand the change of the LDOS in Figs. || (a) to (e). Since the radius 
reduces to in the limit E — > 0, the suppression lines in Fig. || reduce to the lines connecting nearest-neighbor vortex 
centers for E — 0, which is seen in Fig. | (a). 

Another way to examine the quasiparticle excitations in the vortex state is to see how the spectrum evolves along 
the radial directions. We consider the spectrum at the points along the 0° (line OA in Fig. Q), 15° (line OB) 
and 30° (line OC) radial directions. In the figures, we show the spectrum at r = (ca x ,0) for the 0° direction, 
r = (ca x , ca x tan(7r/12)) for the 15° direction and r = {ca xi ca x tan(W6)) for the 30° direction, where we choose the 
point for c=0, 0.1, 0.3, 0.6, 1. As for the case H = 0.1 Tesla, Fig. shows the spectrum at the points along the 0° 
direction. For E < 0.8, the spectrum has the similar structure as that of the single vortex case (see Fig. 2 in Ref. [lj 
for comparison). At higher energy, the peak at E = 1 in the single vortex case becomes broad and shifts a little to 
higher energy side due to the vortex lattice effect. Along the other radial directions, the LDOS has the almost similar 
spectrum as that of the 0° direction. Small differences are shown only for E > 0.8. 

On the other hand, the spectrum for H = 1 Tesla is shown in Fig. [l0| where (a), (b) and (c) are that along 0°, 
15° and 30° radial directions, respectively. It is notable that the lower energy peak in Fig. ^ splits into two or three 
peaks in this higher field case. The behavior of these peaks varies depending on the direction of the radial lines, as 
shown in Fig. [l(] (a)-(c). The peak above E — 1 in Fig. || is suppressed in Fig. [l0| The peak of the spectrum at the 
boundary point (c — 1) shifts to higher energy from 1. 



C. LDOS at higher field 



Here, to discuss the transfer of the quasiparticle bound states between vortex cores, we focus on the LDOS structure 
at E = 0. With increasing a magnetic field, the peak of the LDOS at the vortex center decreases {N{E = 0, r = 0) =71 
for 0.1 Tesla, 65 for 1 Tesla, 44 for 2 Tesla and 11 for 4 Tesla), and the distribution of the LDOS extends sixfold 
symmetrically to the wider region. Therefore, the LDOS around the vortex core is connected each other at high 
magnetic field even in the zero energy state. To show the connection clearly, the LDOS for H — 2 Tesla and 4 Tesla 
is presented in Figs. [□] and [l2| At H = 2 Tesla in Fig. [□] (a), the LDOS around the vortex core is connected each 
other at the middle point of the line AC of Fig. |l|. At further high field H = 4 Tesla shown in Fig. |ll] (b), the LDOS 
is seems to have the component uniformly distributing all over the unit cell. There, the LDOS has small peak at the 
vortex center and minimum at the stationary point of the current flow (point A of Fig. |l|). With raising a magnetic 
field, the LDOS at the boundary of the Wigner-Seitz cell increases (its maximum value is 0.01 for 0.1 Tesla, 0.18 for 1 
Tesla, 0.44 for 2 Tesla and 0.90 for 4 Tesla), which means that the transfer of the quasiparticle bound states between 
vortex cores increases. 

On NbSe2, the dHvA oscillation is observed at magnetic fields down to 4 Tesla in the superconducting mixed state. 
As discussed above, we show that there is large transfer of the bound state between vortex cores in this high field 
region. This transfer seems to be a possible origin of the dHvA oscillation in superconductors. 



IV. SUMMARY AND DISCUSSIONS 



By using the self-consistently obtained pair potential and vector potential, the LDOS of the triangular vortex lattice 
is calculated in an isotropic s-wave superconductors based on the quasi-classical Eilenberger theory. Important results 
obtained for the vortex lattice are as follows; (i) We do find a sixfold star shape of the LDOS and the 30° rotation upon 
elevating the bias energy. The sixfold star originates from the triangular vortex lattice effect, (ii) The orientation 
of the star coincides with the STM data, namely the ray extends toward the next nearest-neighbor vortex direction 
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at lower bias energy. Therefore we succeed in determining the absolute, direction relative to the vortex lattice. This 
is one of the most eminent features in the STM data by Hess et alxfa Thus the experimental features (1) and (2) 
mentioned in Sec. Q are qualitatively reproduced by the naive calculation of the LDOS in the triangular vortex lattice, 
(iii) For E ~ 1, the LDOS distributes around the boundary of the Wigner-Seitz cell of the vortex lattice, and there 
are little distribution around the vortex core region. Therefore, the vortex core is detected as a sixfold symmetric 
dark object in the STM images. At the energy E ~ 1.2, the LDOS has large intensity at the point farthest from the 
vortex center (the point C in Fig. [j]). (iv) The characteristic sixfold symmetric LDOS structure in the low energy 
case appears only at a high magnetic field such as 1 Tesla for the material parameters appropriate to NbSe2, where 
the core regions substantially overlap each other. At a lower magnetic field such as 0.1 Tesla, the LDOS reduces to 
the almost cylindrical structure in the low energy case. 

As mentioned above, the numerically obtained LDOS in the vortex lattice case can qualitatively reproduce the 
characteristic features (1) and (2) in the experimental data on NbSe2- However, the detailed comparison reveals some 
discrepancies between the experimental data and our results, as follows; (v) At the intermediate energy of the 30° 
rotation, our results does not reproduce the split parallel ray structure (the feature (3) mentioned in Sec. |). In 
the experiment, the ray of the star in the 30° direction at E = splits into two parallel rays with elevating E. In 
our results, the ray along the 30° direction at E = weakens and new ray starts to extend along the 0° direction 
with increasing E. (vi) The shape of the star in our results is rather different from the experimental data at E = 0. 
In the experimental data, the ray of the star sharply extends with the narrow width along the 30° direction. On 
the other hand, the ray extends with wide width in our results, (vii) In our results, the structure of the LDOS has 
cylindrical symmetry at a low magnetic field in the low energy case. However, in the experiment by Hess et al., the 
clear star-shaped structure is observed at low fields down to 500 Gauss. 

These discrepancies mean that we have to consider the other effects in addition to the vortex lattice effect studied 
here in order to reproduce the detailed LDOS structure of NbSe2 observed by Hess et al. If we consider the effect of 
the anisotropic s-wave pairing which has sixfold symmetrically anisotropic gap in the basal plane, the STM images 
by Hess et al. are quite well reproduced. We report it in detail elsewhereJla It suggests that the well-studied material 
NbSe2 seems to be one particular system which has anisotropic pairing. Our calculation is performed for the idealized 
isotropic s-wave superconductor. The characteristic features numerically obtained in this paper are expected to be 
well observed in other materials which are isotropic s-wave superconductors. 

To discuss the transfer of the quasiparticle bound states between vortex cores, we consider the LDOS structure at 
E = 0. At high magnetic fields, the LDOS around the vortex core is connected each other even in the zero energy 
state. With raising a magnetic field, this connection increases. It means that the transfer of the quasiparticle bound 
state between vortex cores increases. On NbSe2, the dHvA oscillation is observed at magnetic fields down to 4 Tesla 
in the superconducting mixed state. We show that there is large transfer of the bound state between vortex cores in 
this high field region. This transfer seems to be a possible origin of the dHvA oscillation in superconductors. The 
estimate whether this transfer is enough for the dHvA oscillation or not is a future problem. 
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FIG. 1. The configuration of the vortex lattice in our figures. The vortex centers are shown by •. A hexagon enclosed 
by dashed lines presents the Wigner-Seitz cell of the vortex lattice. When we consider the profile of the spatial variation, we 
present it along the trajectories of line OA (0° radial direction), line OB (15° radial direction), line OC (30° radial direction) 
and line AC (boundary). 

FIG. 2. Spatial variation of the pair potential at T/T c — 0.1. Contour plot of the amplitude, |A(r)|, is presented, (a) 
H = 0.1 Tesla. The region 12. 8£ x 12. 8£ is presented. The core region localizes in a small area around each vortex center, (b) 
H = 1 Tesla. The region 4£ x 4£ is presented, (c) H = 2 Tesla. The region 2.8£ x 2.8£ is presented. For H =1 Tesla and 
2 Tesla, the core region occupies the large part of the unit cell. The pair potential shows cylindrical symmetry at the inner 
region of the vortex core, and sixfold symmetry at the outer region. 

FIG. 3. Spatial variation of the magnetic field at T/T c = 0.1. Contour plot of the normalized value H{r)/H is presented 
for H = 0.1 Tesla (a), 1 Tesla (b) and 2 Tesla (c). For H = 1 Tesla and 2 Tesla, H(r) extends a little to the 30° direction in 
the vortex core region. 

FIG. 4. Spatial variation of the LDOS N(E,r) at H = 0.1 Tesla. The contour plot is shown for £=0.8 (a), 1.0 (b) and 1.2 

(c) . The LDOS structure shows cylindrical symmetry for lower energy E < 0.8, and sixfold symmetry for higher energy E > 1. 

FIG. 5. Spatial variation of the LDOS N(E, r) at H = 0.1 Tesla. The profiles along the lines OA (solid line), OB (dotted 
line), OC (dot-dashed line) and AC (thick line) of Fig. 1 are presented for i?=0.8 (a), 1.0 (b) and 1.2 (c) as a function of r 
(distance from the vortex center) . Each figure corresponds to that of Fig. 4. 

FIG. 6. Spatial variation of the LDOS N(E, r) at H = 1 Tesla. The contour plot is shown for E=0 (a), 0.2 (b), 0.6 (c), 1.0 

(d) and 1.2 (e). In this high magnetic field case, the LDOS shows sixfold symmetry even at the low energy down to E — 0. 



FIG. 7. Spatial variation of the LDOS N(E,r) at H = 1 Tesla. The profiles along the lines OA (solid line), OB (dotted 
line), OC (dot-dashed line) and AC (thick line) of Fig. 1 are presented for E=0 (a), 0.2 (b), 0.6 (c), 1.0 (d) and 1.2 (e) as a 
function of r. Each figure corresponds to that of Fig. 6. 

FIG. 8. Schematically presented LDOS structure, which corresponds to that of Fig. 6. Without the vortex lattice effect, 
the LDOS has the cylindrical structure, which is presented as white rings. To this structure, the vortex lattice effect suppresses 
the LDOS along the lines presented in the figure. Therefore, it is seen that the LDOS in Fig. 6 has small value at the points 
along the lines. 

FIG. 9. Spectrum N(E,r) at the points r along the 0° radial direction line OA for H = 0.1 Tesla. The sample points are 
given by r = (ca x ,0) with c=0, 0.1, 0.3, 0.6 and 1. 

FIG. 10. Spectrum N(E,r) at the points r along the 0° radial direction line OA (a), the 15° line OB (b) and the 30° line 
OC (c) for H — 1 Tesla. The sample points are given by r = (ca x , 0) for (a), (ca x ,ca x tan(7r/12)) for (b) and (ca x , ca x tan(7r/6)) 
for (c), with c=0, 0.1, 0.3, 0.6 and 1. 

FIG. 11. Spatial variation of the LDOS N(E, r) for E = 0. Contour plot is shown for H =2 Tesla (a) and 4 Tesla (b). At 
higher magnetic fields, the LDOS around the vortex is connected each other. 

FIG. 12. Spatial variation of the LDOS N(E,r) for E = 0. Profile along the line OA (solid line), OB (dotted line), OC 
(dot-dashed line) and AC (thick line) of Fig. 1. are presented for H —2 Tesla (a) and 4 Tesla (b) as a function of r. Each 
figure corresponds to that of Fig. 11. 
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